In [2]:
using Random, DelimitedFiles, PyPlot, StatsBase

Load reference data and define plotting (this will change after Przemek prepares better visualization)

In [3]:
loc = readdlm("loc.txt")
Out[3]:
153×2 Array{Float64,2}:
 -79.2909  43.7837
 -79.0132  43.8662
 -79.8246  44.1615
 -83.109   42.0969
 -80.0141  43.1983
 -79.4112  44.015 
 -75.742   45.2714
 -79.6861  44.3328
 -79.7098  44.4104
 -77.3964  44.1944
 -79.7113  43.864 
 -78.7147  43.9102
 -79.2925  45.0493
   â‹®              
 -79.5551  43.8455
 -82.385   42.5862
 -80.0065  44.5249
 -79.9096  43.3162
 -80.5485  43.509 
 -80.5151  43.4696
 -79.2294  43.0148
 -78.9597  43.9174
 -83.011   42.2707
 -82.9424  42.314 
 -79.6284  43.7861
 -80.7346  43.1151
In [4]:
dist = readdlm("distkm.txt")
Out[4]:
153×153 Array{Float64,2}:
   0.0      25.5691   81.8831  388.798   …  366.945    32.1738  147.782 
  25.1537    0.0     102.777   413.394      391.541    56.0079  172.378 
  82.3007  103.224     0.0     399.893      378.04     50.7209  158.877 
 388.747   413.54    399.865     0.0         36.9392  367.115   241.551 
  96.4183  117.187   129.688   303.684      281.831    82.3623   62.6679
  31.5325   50.8019   52.7316  406.836   …  384.982    39.8438  165.82  
 359.209   334.3     401.386   745.918      724.065   381.755   504.902 
  83.9051  104.139    29.3792  428.76       406.907    65.2319  187.744 
  89.3415  109.575    36.8667  436.247      414.394    72.7194  195.231 
 163.578   138.978   238.802   551.408      529.554   194.986   310.392 
  44.8736   67.9438   38.4316  371.541   …  349.688    13.2937  130.525 
  50.7736   25.6199  125.444   438.685      416.832    81.6278  197.669 
 166.597   149.768   124.579   523.959      502.106   160.431   282.943 
   ⋮                                     ⋱    ⋮                         
  30.3718   52.6      52.9476  378.209      356.355    11.1236  137.193 
 308.335   333.129   318.098   118.835       93.3661  286.704   162.173 
 122.688   142.922    50.7985  435.24       413.008   100.685   194.224 
  78.0129   99.4273  111.722   317.18       295.327    63.957    76.1638
 121.065   145.858   109.244   300.795   …  277.808    99.0817   61.2077
 116.377   141.17    112.958   295.812      272.824    94.7451   56.111 
 143.864   163.623   178.233   366.067      344.213   129.714   136.748 
  34.8408   10.1703  103.898   421.578      399.724    59.6691  180.562 
 370.41    395.203   381.528    26.0345      12.0937  348.778   223.214 
 366.305   391.098   377.423    35.396   …    0.0     344.673   219.109 
  32.2039   56.4238   50.6923  367.576      345.723     0.0     126.56  
 148.379   173.172   159.497   241.625      219.771   126.747     0.0   
In [5]:
function plot_tour(t, loc)
    t2 = [t; t[1]]
    plot(loc[t2, 1], loc[t2, 2])
end
Out[5]:
plot_tour (generic function with 1 method)
In [6]:
#This code is using route data files that need to be downloaded and unzipped from:
#https://github.com/pszufe/MapDatasets/raw/master/Ontario/walmart_routing.zip

pos = Dict{Int,Tuple{Float64,Float64}}()
for line in readlines("walmart_node_locations.txt")
    els = split(line,"\t")
    pos[parse(Int,els[1])] = (parse(Float64,els[2]),parse(Float64,els[3]))
end

routes = Matrix{Vector{Int}}(undef,size(dist)...)
fill!(routes,Int[])
f = open("walmart_routes.txt")
while !eof(f)
    line = split(readline(f),"\t")
    i,j = parse.(Int,line[2:3])
    routes[i,j]= parse.(Int,split(line[12],"#"))
end
close(f)
In [7]:
k = first(keys(pos))
const MAP_BOUNDS = [[pos[k]...],[pos[k]...]]
for v in values(pos)
    (v[1] < MAP_BOUNDS[1][1]) && (MAP_BOUNDS[1][1]=v[1])
    (v[2] < MAP_BOUNDS[1][2]) && (MAP_BOUNDS[1][2]=v[2])
    (v[1] > MAP_BOUNDS[2][1]) && (MAP_BOUNDS[2][1]=v[1])
    (v[2] > MAP_BOUNDS[2][2]) && (MAP_BOUNDS[2][2]=v[2])
    
end
MAP_BOUNDS
Out[7]:
2-element Array{Array{Float64,1},1}:
 [42.0339, -94.4713]
 [49.8608, -74.6007]
In [8]:
#required installation for map vizualiztion
#using Conda
#Conda.runconda(`install folium -c conda-forge`)
using PyCall
flm = pyimport("folium")

function plot_tour2(t)
    m = flm.Map()
    for k=1:length(t)
        i = t[k]
        j = t[k<length(t) ? (k+1) : 1]
        info = "Route $(k)-th <br>from=$i to=$j\n<br>distance=$(round(dist[i,j],digits=3))"
        flm.PolyLine(
            [pos[n] for n in routes[i,j]],
            popup=info,
            tooltip=info
        ).add_to(m)
        
    end
       

    for k=1:length(t)
        i = t[k]
        j = t[k<length(t) ? (k+1) : 1]
        info = "Node $(k)-th <br>number=$i next=$j\n<br>\n$(round.(pos[routes[i,j][1]],digits=6))"
        flm.Circle(
          location=pos[routes[i,j][1]],
          popup=info,
          tooltip=info,
          radius=2000,
          color="crimson",
          fill=true,
          fill_color="crimson"
       ).add_to(m)
    end
    
    m.fit_bounds(Tuple.(MAP_BOUNDS))
    m
end
Out[8]:
plot_tour2 (generic function with 1 method)

Tour length calculation function common for all algorithms

In [9]:
function eval_tour(tour, d)
    travel = d[tour[end], tour[1]]
    for i in 2:length(tour)
        travel += d[tour[i-1], tour[i]]
    end
    travel
end
Out[9]:
eval_tour (generic function with 1 method)

Random initial tour generation common for all algorithms

In [10]:
init_tour(n) = randperm(n)
Out[10]:
init_tour (generic function with 1 method)

Function choosing one random neighbor used in local search and SANN

In [11]:
function pick_neighbor(tour, prob)
    neighbor = copy(tour)

    x, y = rand(eachindex(tour), 2)
    if x != y
        if rand() < prob
            x,y = minmax(x,y)
            reverse!(@view neighbor[x:y])
        else
            tmp = neighbor[x]
            if x < y
                neighbor[x:y-1] = @view neighbor[x+1:y]
            else
                neighbor[y+1:x] = @view neighbor[y:x-1]
            end
            neighbor[y] = tmp
        end
    end
    return neighbor
end
Out[11]:
pick_neighbor (generic function with 1 method)

Local search routine

In [12]:
function localsearch(dist, iter, prob)
    t = init_tour(size(dist, 1))
    v = eval_tour(t, dist)
    best_t, best_v = t, v
    for i in 1:iter
        newt = pick_neighbor(t, prob)
        newv = eval_tour(newt, dist)
        if newv < v
            t, v = newt, newv
            if v < best_v
                best_v, best_t = v, t
            end
        end
    end
    (best_v, best_t)
end
Out[12]:
localsearch (generic function with 1 method)
In [13]:
alg_localsearch(;iter=10^7, prob=0.8) = (x -> localsearch(x, iter, prob), "Local search (iter=$iter, prob=$prob)")
Out[13]:
alg_localsearch (generic function with 1 method)

Simulated annealing routine

In [14]:
function sann(dist, T0, Tf, α, iter, prob)
    T = T0
    t = init_tour(size(dist, 1))
    v = eval_tour(t, dist)
    best_t, best_v = t, v
    while T >= Tf
        for i in 1:iter
            newt = pick_neighbor(t, prob)
            newv = eval_tour(newt, dist)
            if newv < v || rand() < exp(-(newv - v) / T)
                t, v = newt, newv
                if v < best_v
                    best_v, best_t = v, t
                end
            end
        end
        T *= α
    end
    (best_v, best_t)
end
Out[14]:
sann (generic function with 1 method)
In [15]:
alg_sann(; T0=1, Tf=0.001, α=0.9, iter=round(10^7 * log(α)/ log(Tf/T0)), prob=0.8) =
    (x -> sann(x, T0, Tf, α, iter, prob), "SANN (T0=$T0, Tf=$Tf, α=$α, iter=$iter, prob=$prob)")
Out[15]:
alg_sann (generic function with 1 method)

Genetic algorithm routine

In [16]:
function init_tour_ga(d)
    t = randperm(size(d, 1))
    (eval_tour(t, d), t)
end
Out[16]:
init_tour_ga (generic function with 1 method)
In [17]:
function mutate(tour, d)
    x, y = minmax(rand(eachindex(tour[2]), 2)...)
    reverse!(@view tour[2][x:y])
    return (eval_tour(tour[2], d), tour[2])
end
Out[17]:
mutate (generic function with 1 method)
In [18]:
function crossover(t1, t2, d)
    x, y = minmax(rand(eachindex(t1[2]), 2)...)
    child = similar(t1[2])
    used = Set{Int}()
    for i in x:y
        child[i] = t1[2][i]
        push!(used, t1[2][i])
    end
    pos = 1
    for j in eachindex(t2[2])
        if pos == x
            pos = y + 1
        end
        if !(t2[2][j] in used)
            child[pos] = t2[2][j]
            pos += 1
        end
    end
    (eval_tour(child, d), child)
end
Out[18]:
crossover (generic function with 1 method)
In [19]:
function ga(dist, popsize, iter, mutprob)
    w_replace = Weights(0:popsize-1)
    w_parent = Weights(popsize:-1:1)
    pop = [init_tour_ga(dist) for i in 1:popsize]
    sort!(pop)
    best_v, best_t = pop[1]
    for i in 1:iter
        replace_idx = sample(1:popsize, w_replace)
        pop[replace_idx] = crossover(pop[sample(1:popsize, w_parent)],
                                     pop[sample(1:popsize, w_parent)],
                                     dist)
        for j in 2:popsize
            if rand() < mutprob
                pop[j] = mutate(pop[j], dist)
            end
        end
        sort!(pop)
        if pop[1][1] < best_v
            best_v, best_t = pop[1]
        end
    end
    (best_v, best_t)
end
Out[19]:
ga (generic function with 1 method)
In [20]:
alg_ga(; popsize=120, iter=3*10^5, mutprob=0.01) =
    (x -> ga(x, popsize, iter, mutprob), "GA popsize=$popsize, iter=$iter, mutprob=$mutprob")
Out[20]:
alg_ga (generic function with 1 method)

Tabu search routine

In [21]:
function list_neighbors(tour, d)
    neighbors = Vector{typeof(tour)}()
    n = length(tour)

    for x in 1:n-1, y in x+1:n
        neighbor = copy(tour)
        reverse!(@view neighbor[x:y])
        push!(neighbors, neighbor)
    end
    sort!([(eval_tour(nei, d), nei) for nei in neighbors])
end
Out[21]:
list_neighbors (generic function with 1 method)
In [22]:
function tabusearch(d, tabulength, iter)
    t = init_tour(size(dist, 1))
    v = eval_tour(t, d)
    best_t, best_v = t, v
    tabulist = [Int[] for i in 1:tabulength]
    tabupos = mod1(2, tabulength)
    tabulist[tabupos] = t
    for i in 1:iter
        neighbors = list_neighbors(t, d)
        found = false
        for (neiv, nei) in neighbors
            if !(nei in tabulist)
                found = true
                tabulist[tabupos] = nei
                tabupos = mod1(tabupos+1, tabulength)
                t, v = nei, neiv
                if v < best_v
                    best_v, best_t = v, t
                end
                break
            end
        end
        if !found
            break
        end
    end
    (best_v, best_t)
end
Out[22]:
tabusearch (generic function with 1 method)
In [23]:
alg_tabusearch(; tabulength=100, iter=1000) =
    (x -> tabusearch(x, tabulength, iter), "Tabu search tabulength=$tabulength iter=$iter")
Out[23]:
alg_tabusearch (generic function with 1 method)

Algorithm comparison

In [24]:
function run(;loc, dist, alg)
    @time res = alg[1](dist)
    #plot_tour(res[2], loc)
    #title("$(alg[2]); cost: $(res[1])")
    plot_tour2(res[2])    
end
Out[24]:
run (generic function with 1 method)
In [25]:
Random.seed!(2019);
In [26]:
run(loc=loc, dist=dist, alg=alg_localsearch())
  8.079232 seconds (33.91 M allocations: 15.056 GiB, 11.15% gc time)
Out[26]:
┌ Warning: `getindex(o::PyObject, s::AbstractString)` is deprecated in favor of dot overloading (`getproperty`) so elements should now be accessed as e.g. `o."s"` instead of `o["s"]`.
│   caller = show(::IOContext{Base.GenericIOBuffer{Array{UInt8,1}}}, ::MIME{Symbol("text/html")}, ::PyObject) at PyCall.jl:888
â”” @ PyCall C:\JuliaPkg\Julia-1.0.3\packages\PyCall\RQjD7\src\PyCall.jl:888
In [27]:
run(loc=loc, dist=dist, alg=alg_sann())
  7.899815 seconds (34.13 M allocations: 15.156 GiB, 12.31% gc time)
Out[27]:
In [28]:
run(loc=loc, dist=dist, alg=alg_ga())
  4.338142 seconds (5.70 M allocations: 1.410 GiB, 4.10% gc time)
Out[28]:
In [29]:
run(loc=loc, dist=dist, alg=alg_tabusearch())
 12.664319 seconds (34.90 M allocations: 15.969 GiB, 21.21% gc time)
Out[29]: